#include "MeshCutter.h"
#include <chrono>
#include <unordered_set>
#include <algorithm>
#include <Eigen/Dense>
#include <sstream>

void MeshCutter::processSurface(const std::vector<std::shared_ptr<ImplicitSurface>>& surfaces, 
                               int surfaceType) {
    for (const auto& surface : surfaces) {
        processSurface(*surface, surfaceType);
    }
}

void MeshCutter::processSurface(const ImplicitSurface& surface, int surfaceType) {
    // 如果是第一个地层界面(11)，将所有z坐标最小的点标记为该界面点
    if (surfaceType == 11) {
        double minZ = nodes_[0].z;
        // 先找到最小的z坐标
        for (const auto& node : nodes_) {
            minZ = std::min(minZ, node.z);
        }
        // 将所有z坐标等于minZ的点标记为界面点
        for (auto& node : nodes_) {
            if (std::abs(node.z - minZ) < 1e-10) {
                node.type = surfaceType;
            }
        }
        return;
    }
    
    // 如果是最后一个地层界面，将所有z坐标最大的点标记为该界面点
    if (surfaceType == lastStrataType_) {
        double maxZ = nodes_[0].z;
        // 先找到最大的z坐标
        for (const auto& node : nodes_) {
            maxZ = std::max(maxZ, node.z);
        }
        // 将所有z坐标等于maxZ的点标记为界面点
        for (auto& node : nodes_) {
            if (std::abs(node.z - maxZ) < 1e-10) {
                node.type = surfaceType;
            }
        }
        return;
    }
    
    // 添加预处理步骤
    preprocessStrataInterface(surface, surfaceType);
    
    auto start = std::chrono::high_resolution_clock::now();
    
    // 获取界面的包围盒
    const auto& bounds = surface.getBounds();
    
    // 找到包围盒内的所有边
    std::unordered_set<std::pair<size_t, size_t>, PairHash> edgesInBounds;
    for (const auto& tet : elements_) {
        std::vector<std::pair<size_t, size_t>> edges = {
            {tet.v1, tet.v2}, {tet.v1, tet.v3}, {tet.v1, tet.v4},
            {tet.v2, tet.v3}, {tet.v2, tet.v4}, {tet.v3, tet.v4}
        };
        
        for (auto& edge : edges) {
            // 标准化边的顶点顺序
            if (edge.first > edge.second) {
                std::swap(edge.first, edge.second);
            }
            
            const auto& p1 = nodes_[edge.first];
            const auto& p2 = nodes_[edge.second];
            
            // 检查边是否在包围盒内
            if (isEdgeInBounds(p1, p2, bounds.first, bounds.second)) {
                edgesInBounds.insert(edge);
            }
        }
    }
    
    // 找到需要切割的边
    std::vector<std::pair<size_t, size_t>> activeEdges;
    for (const auto& edge : edgesInBounds) {
        const auto& p1 = nodes_[edge.first];
        const auto& p2 = nodes_[edge.second];
        
        // 如果任一端点已经是界面点，跳过
        if (p1.type == surfaceType || p2.type == surfaceType) {
            continue;
        }
        
        double phi1 = surface.evaluatePhi(p1);
        double phi2 = surface.evaluatePhi(p2);
        
        // 只有当两个端点都不在界面上，且在界面的不同侧时，才需要切割
        if (std::abs(phi1) >= 1e-10 && std::abs(phi2) >= 1e-10 && phi1 * phi2 < 0) {
            activeEdges.push_back(edge);
        }
    }
    
    // 如果没有需要切割的边，直接返回
    if (activeEdges.empty()) {
        return;
    }
    
    // 计算切割点并切割网格
    auto cutPoints = computeCutPoints(activeEdges, surface, surfaceType);
    
    // 更新统计信息
    stats_[surfaceType].cutEdges = static_cast<int>(activeEdges.size());
    stats_[surfaceType].newNodes = static_cast<int>(cutPoints.size());
    
    // 切割四面体
    splitTetrahedra(cutPoints, surfaceType);
    
    // 更新网格数据结构，为下一个界面做准备
    updateMeshData();
    
    // 计算处理时间
    auto end = std::chrono::high_resolution_clock::now();
    stats_[surfaceType].time = std::chrono::duration<float>(end - start).count();
}

std::vector<std::pair<size_t, size_t>> MeshCutter::findActiveEdges(
    const ImplicitSurface& surface,
    int surfaceType) {
    
    std::vector<std::pair<size_t, size_t>> activeEdges;
    
    // 标记所有在界面上的点
    for (size_t i = 0; i < nodes_.size(); ++i) {
        if (std::abs(surface.evaluatePhi(nodes_[i])) < 1e-10) {
            nodes_[i].type = surfaceType;  // 直接修改原节点类型
        }
    }
    
    // 获取界面的包围盒
    const auto& bounds = surface.getBounds();
    
    // 找到需要切割的边
    for (const auto& tet : elements_) {
        std::vector<std::pair<size_t, size_t>> edges = {
            {tet.v1, tet.v2}, {tet.v1, tet.v3}, {tet.v1, tet.v4},
            {tet.v2, tet.v3}, {tet.v2, tet.v4}, {tet.v3, tet.v4}
        };
        
        for (auto& edge : edges) {
            // 标准化边的顶点顺序
            if (edge.first > edge.second) {
                std::swap(edge.first, edge.second);
            }
            
            const auto& p1 = nodes_[edge.first];
            const auto& p2 = nodes_[edge.second];
            
            // 如果任一端点已经是界面点，跳过
            if (p1.type == surfaceType || p2.type == surfaceType) {
                continue;
            }
            
            // 检查边是否在包围盒内
            if (!isEdgeInBounds(p1, p2, bounds.first, bounds.second)) {
                continue;
            }
            
            // 对于断层界面，检查边界函数
            if (surfaceType >= 1 && surfaceType <= 10) {  // 断层界面
                if (surface.evaluatePsi(p1) > 0 && surface.evaluatePsi(p2) > 0) {
                    continue;  // 如果两个端点都在断层边界外，跳过这条边
                }
            }
            
            double phi1 = surface.evaluatePhi(p1);
            double phi2 = surface.evaluatePhi(p2);
            
            // 只有当两个端点都不在界面上，且在界面的不同侧时，才需要切割
            if (std::abs(phi1) >= 1e-10 && std::abs(phi2) >= 1e-10 && phi1 * phi2 < 0) {
                activeEdges.push_back(edge);
            }
        }
    }
    
    return activeEdges;
}

std::unordered_map<std::pair<size_t, size_t>, size_t, PairHash>
MeshCutter::computeCutPoints(
    const std::vector<std::pair<size_t, size_t>>& activeEdges,
    const ImplicitSurface& surface,
    int surfaceType) {
    
    std::unordered_map<std::pair<size_t, size_t>, size_t, PairHash> cutPoints;
    
    for (const auto& edge : activeEdges) {
        const auto& p1 = nodes_[edge.first];
        const auto& p2 = nodes_[edge.second];
        
        // 计算切割点
        Point3D cutPoint = computeCutPoint(p1, p2, surface);
        
        // 对于断层界面，检查边界
        if (surfaceType >= 1 && surfaceType <= 10) {  // 断层界面
            double psi = surface.evaluatePsi(cutPoint);
            if (psi > 0) {
                continue;  // 如果切割点在断层边界外，跳过
            }
        }
        
        // 添加新的界面点
        cutPoint.type = surfaceType;
        nodes_.push_back(cutPoint);
        cutPoints[edge] = nodes_.size() - 1;
    }
    
    return cutPoints;
}

void MeshCutter::splitTetrahedra(
    const std::unordered_map<std::pair<size_t, size_t>, size_t, PairHash>& cutPoints,
    int surfaceType) {
    
    std::vector<Tetrahedron> newElements;
    newElements.reserve(elements_.size() * 2);
    
    for (const auto& tet : elements_) {
        // 获取当前四面体的切割点
        std::unordered_map<Edge, size_t> tetCuts;
        std::vector<Edge> tetEdges = {
            Edge(tet.v1, tet.v2), Edge(tet.v1, tet.v3),
            Edge(tet.v1, tet.v4), Edge(tet.v2, tet.v3),
            Edge(tet.v2, tet.v4), Edge(tet.v3, tet.v4)
        };
        
        // 收集切割点
        for (const auto& edge : tetEdges) {
            std::pair<size_t, size_t> key(edge.v1, edge.v2);
            if (key.first > key.second) std::swap(key.first, key.second);
            
            auto it = cutPoints.find(key);
            if (it != cutPoints.end()) {
                tetCuts[edge] = it->second;
            }
        }
        
        // 根据切割边数量选择切割方案
        std::vector<Tetrahedron> splitResult;
        switch (tetCuts.size()) {
            case 4: splitResult = splitCaseA(tet, tetCuts); break;
            case 3: splitResult = splitCaseB(tet, tetCuts); break;
            case 2: splitResult = splitCaseC(tet, tetCuts); break;
            case 1: splitResult = splitCaseD(tet, tetCuts); break;
            default: splitResult = {tet}; break;
        }
        
        newElements.insert(newElements.end(), splitResult.begin(), splitResult.end());
    }
    
    elements_ = std::move(newElements);
    updateMeshData();
}

bool MeshCutter::isPointOnSurface(const Point3D& p, const ImplicitSurface& surface) const {
    return std::abs(surface.evaluatePhi(p)) < 1e-10;
}

Point3D MeshCutter::computeCutPoint(
    const Point3D& p1, const Point3D& p2,
    const ImplicitSurface& surface) const {
    
    double phi1 = surface.evaluatePhi(p1);
    double phi2 = surface.evaluatePhi(p2);
    
    // 使用线性插值计算切割点
    double t = std::abs(phi1) / (std::abs(phi1) + std::abs(phi2));
    
    return Point3D(
        p1.x + t * (p2.x - p1.x),
        p1.y + t * (p2.y - p1.y),
        p1.z + t * (p2.z - p1.z)
    );
}

int MeshCutter::computeCutPattern(
    const Tetrahedron& tet,
    const std::unordered_map<std::pair<size_t, size_t>, size_t, PairHash>& cutPoints) const {
    
    int pattern = 0;
    
    // 按照Python代码的顺序定义边
    // 注意：这里的顺序必须和Python代码完全一致
    std::vector<std::pair<size_t, size_t>> edges = {
        {tet.v1, tet.v2}, {tet.v1, tet.v3}, {tet.v1, tet.v4},  // v1的边
        {tet.v2, tet.v3}, {tet.v2, tet.v4},                     // v2的边
        {tet.v3, tet.v4}                                        // v3的边
    };
    
    // 检查每条边是否被切割
    for (size_t i = 0; i < edges.size(); ++i) {
        auto edge = edges[i];
        // 标准化边的顶点顺序
        if (edge.first > edge.second) {
            std::swap(edge.first, edge.second);
        }
        
        // 检查边是否在切割点集合中
        if (cutPoints.find(edge) != cutPoints.end()) {
            pattern |= (1 << i);
        }
    }
    
    return pattern;
}

std::vector<Tetrahedron> MeshCutter::generateNewTetrahedra(
    const Tetrahedron& tet,
    const std::unordered_map<std::pair<size_t, size_t>, size_t, PairHash>& cutPoints,
    int pattern) const {
    
    // 将pair类型的cutPoints转换为Edge类型
    std::unordered_map<Edge, size_t> edgeCuts;
    for (const auto& [edge, idx] : cutPoints) {
        // 标准化边的顶点顺序
        size_t v1 = edge.first, v2 = edge.second;
        if (v1 > v2) std::swap(v1, v2);
        Edge e(v1, v2);
        
        // 检查是否是当前四面体的边
        if (pattern & (1 << getEdgeIndex(e))) {
            edgeCuts[e] = idx;
        }
    }
    
    // 获取四面体的顶点（按照Python代码的顺序）
    // Python代码中：a,b,c,d = tet[3],tet[0],tet[1],tet[2]
    size_t a = tet.v4;  // 第四个顶点
    size_t b = tet.v1;  // 第一个顶点
    size_t c = tet.v2;  // 第二个顶点
    size_t d = tet.v3;  // 第三个顶点
    
    // 获取所有可能的边
    std::unordered_map<std::string, Edge> edges = {
        {"ab", Edge(a, b)}, {"ac", Edge(a, c)}, {"ad", Edge(a, d)},
        {"bc", Edge(b, c)}, {"bd", Edge(b, d)}, {"cd", Edge(c, d)}
    };
    
    // 检查哪些边被切割
    std::unordered_set<std::string> cutEdges;
    for (const auto& [name, edge] : edges) {
        if (edgeCuts.count(edge)) {
            cutEdges.insert(name);
        }
    }
    
    // 根据切割边的组合选择切割方案
    if (cutEdges.size() == 4) {
        // 严格检查每种模式的边组合
        const std::vector<std::vector<std::string>> patterns = {
            {"ab", "ac", "cd", "bd"},  // A1
            {"ad", "bd", "bc", "ac"},  // A2
            {"ab", "bc", "cd", "ad"}   // A3
        };
        
        for (const auto& pattern : patterns) {
            if (cutEdges == std::unordered_set<std::string>(pattern.begin(), pattern.end())) {
                return splitCaseA(tet, edgeCuts);
            }
        }
    }
    else if (cutEdges.size() == 3) {
        // Case B: 三条边被切割
        std::vector<std::string> patterns[] = {
            {"ab", "ac", "ad"},  // B1
            {"ad", "cd", "bd"},  // B2
            {"ab", "bd", "bc"},  // B3
            {"ac", "bc", "cd"}   // B4
        };
        
        for (const auto& pattern : patterns) {
            bool match = true;
            for (const auto& edge : pattern) {
                if (!cutEdges.count(edge)) {
                    match = false;
                    break;
                }
            }
            if (match) {
                return splitCaseB(tet, edgeCuts);
            }
        }
    }
    else if (cutEdges.size() == 2) {
        if (cutEdges.count("ab") && cutEdges.count("ac")) {
            return splitCaseC(tet, edgeCuts);  // Case C1
        }
        if (cutEdges.count("bc") && cutEdges.count("cd")) {
            return splitCaseC(tet, edgeCuts);  // Case C2
        }
        if (cutEdges.count("ac") && cutEdges.count("ad")) {
            return splitCaseC(tet, edgeCuts);  // Case C3
        }
        if (cutEdges.count("ab") && cutEdges.count("ad")) {
            return splitCaseC(tet, edgeCuts);  // Case C4
        }
    }
    else if (cutEdges.size() == 1) {
        for (const auto& edge : {"ab", "cd", "bd", "bc", "ad", "ac"}) {
            if (cutEdges.count(edge)) {
                return splitCaseD(tet, edgeCuts);  // Case D
            }
        }
    }
    
    return {tet};  // 如果没有匹配的切割模式，返回原四面体
}

// 辅助函数：获取边的索引
int MeshCutter::getEdgeIndex(const Edge& edge) const {
    // 标准化边的顶点顺序
    size_t v1 = edge.v1, v2 = edge.v2;
    if (v1 > v2) std::swap(v1, v2);
    
    // 直接使用边的顶点索引来确定边的位置
    // 按照Python代码的边索引顺序：[v1-v2, v1-v3, v1-v4, v2-v3, v2-v4, v3-v4]
    std::vector<std::pair<size_t, size_t>> standardEdges = {
        {1, 2}, {1, 3}, {1, 4},  // v1的边
        {2, 3}, {2, 4},          // v2的边
        {3, 4}                   // v3的边
    };
    
    // 查找边在标准顺序中的位置
    for (size_t i = 0; i < standardEdges.size(); ++i) {
        if (v1 == standardEdges[i].first && v2 == standardEdges[i].second) {
            return static_cast<int>(i);
        }
    }
    
    return -1;  // 不应该发生
}

std::unordered_set<Edge> MeshCutter::buildEdgesForTet(const Tetrahedron& tet) const {
    std::unordered_set<Edge> edges;
    edges.insert(Edge(tet.v1, tet.v2));
    edges.insert(Edge(tet.v1, tet.v3));
    edges.insert(Edge(tet.v1, tet.v4));
    edges.insert(Edge(tet.v2, tet.v3));
    edges.insert(Edge(tet.v2, tet.v4));
    edges.insert(Edge(tet.v3, tet.v4));
    return edges;
}

std::vector<size_t> MeshCutter::findTetsForEdge(const Edge& edge) const {
    std::vector<size_t> tetIndices;
    for (size_t i = 0; i < elements_.size(); ++i) {
        const auto& tet = elements_[i];
        auto edges = buildEdgesForTet(tet);
        if (edges.count(edge)) {
            tetIndices.push_back(i);
        }
    }
    return tetIndices;
}

void MeshCutter::splitTets(const std::unordered_map<Edge, size_t>& cutPoints, int surfaceType) {
    const size_t batchSize = 500;  // 使用更小的批处理大小
    
    std::vector<Tetrahedron> newTets;
    std::unordered_set<size_t> processedTets;
    
    try {
        // 分批收集需要处理的四面体
        std::vector<size_t> affectedTets;
        affectedTets.reserve(1000);  // 使用固定的小预分配
        
        for (const auto& [edge, _] : cutPoints) {
            const auto& tetIndices = edgeToTets_[edge];
            affectedTets.insert(affectedTets.end(), tetIndices.begin(), tetIndices.end());
            
            // 如果大小超过阈值，进行一次去重
            if (affectedTets.size() > 5000) {
                std::sort(affectedTets.begin(), affectedTets.end());
                affectedTets.erase(std::unique(affectedTets.begin(), affectedTets.end()), affectedTets.end());
            }
        }
        
        // 最终去重
        std::sort(affectedTets.begin(), affectedTets.end());
        affectedTets.erase(std::unique(affectedTets.begin(), affectedTets.end()), affectedTets.end());
        
        // 使用固定大小的临时缓冲区
        std::vector<Tetrahedron> tempTets;
        tempTets.reserve(batchSize * 2);
        
        // 分批处理
        for (size_t i = 0; i < affectedTets.size(); i += batchSize) {
            tempTets.clear();
            size_t end = std::min(i + batchSize, affectedTets.size());
            
            for (size_t j = i; j < end; ++j) {
                size_t tetIdx = affectedTets[j];
                if (processedTets.count(tetIdx)) continue;
                
                const auto& tet = elements_[tetIdx];
                
                // 获取四面体的切割点
                std::unordered_map<Edge, size_t> tetCuts;
                std::vector<Edge> tetEdges = {
                    Edge(tet.v1, tet.v2), Edge(tet.v1, tet.v3),
                    Edge(tet.v1, tet.v4), Edge(tet.v2, tet.v3),
                    Edge(tet.v2, tet.v4), Edge(tet.v3, tet.v4)
                };
                
                for (const auto& edge : tetEdges) {
                    auto it = cutPoints.find(edge);
                    if (it != cutPoints.end()) {
                        tetCuts[edge] = it->second;
                    }
                }
                
                // 处理四面体
                std::vector<Tetrahedron> splitResult;
                switch (tetCuts.size()) {
                    case 4: splitResult = splitCaseA(tet, tetCuts); break;
                    case 3: splitResult = splitCaseB(tet, tetCuts); break;
                    case 2: splitResult = splitCaseC(tet, tetCuts); break;
                    case 1: splitResult = splitCaseD(tet, tetCuts); break;
                    default: splitResult = {tet}; break;
                }
                
                tempTets.insert(tempTets.end(), splitResult.begin(), splitResult.end());
                processedTets.insert(tetIdx);
            }
            
            // 合并结果
            if (!tempTets.empty()) {
                if (newTets.empty()) {
                    newTets = std::move(tempTets);
                } else {
                    newTets.insert(newTets.end(), tempTets.begin(), tempTets.end());
                }
            }
        }
        
        // 更新网格
        std::vector<Tetrahedron> updatedTets;
        updatedTets.reserve(elements_.size() - processedTets.size() + newTets.size());
        
        for (size_t i = 0; i < elements_.size(); ++i) {
            if (!processedTets.count(i)) {
                updatedTets.push_back(elements_[i]);
            }
        }
        updatedTets.insert(updatedTets.end(), newTets.begin(), newTets.end());
        
        elements_ = std::move(updatedTets);
        stats_[surfaceType].newTets = static_cast<int>(newTets.size() - elements_.size());
        
    } catch (const std::bad_alloc&) {
        // 如果内存分配失败，保持原始网格不变
    }
}

std::vector<Tetrahedron> MeshCutter::splitCaseA(
    const Tetrahedron& tet,
    const std::unordered_map<Edge, size_t>& cuts) const {
    
    try {
        // 获取原始顶点
        size_t b = tet.v1, c = tet.v2, d = tet.v3, a = tet.v4;
        
        // 获取所有可能的边
        std::unordered_map<std::string, Edge> edges = {
            {"ab", Edge(a, b)}, {"ac", Edge(a, c)}, {"ad", Edge(a, d)},
            {"bc", Edge(b, c)}, {"bd", Edge(b, d)}, {"cd", Edge(c, d)}
        };
        
        // 检查哪些边被切割
        std::unordered_set<std::string> cutEdges;
        for (const auto& [name, edge] : edges) {
            if (cuts.count(edge)) {
                cutEdges.insert(name);
            }
        }
        
        // 情况1：ab,ac,cd,bd被切割
        if (cutEdges == std::unordered_set<std::string>{"ab","ac","cd","bd"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            size_t f = cuts.at(edges["ac"]);  // ac的切割点
            size_t g = cuts.at(edges["cd"]);  // cd的切割点
            size_t h = cuts.at(edges["bd"]);  // bd的切割点
            
            // 检查切割点是否共面
            std::vector<Point3D> points = {nodes_[e], nodes_[f], nodes_[g], nodes_[h]};
            if (!checkCoplanar(points)) {
                return {tet};
            }
            
            return {
                // 三棱柱aef-dhg: acb-d方式切分
                Tetrahedron{a, f, e, d},
                Tetrahedron{d, h, f, e},
                Tetrahedron{d, h, g, f},
                // 三棱柱gfc-heb: acb-d方式切分
                Tetrahedron{g, c, f, h},
                Tetrahedron{h, b, c, f},
                Tetrahedron{h, e, b, f}
            };
        }
        
        // 情况2：ad,bd,bc,ac被切割
        else if (cutEdges == std::unordered_set<std::string>{"ad","bd","bc","ac"}) {
            size_t e = cuts.at(edges["ad"]);  // ad的切割点
            size_t f = cuts.at(edges["bd"]);  // bd的切割点
            size_t g = cuts.at(edges["bc"]);  // bc的切割点
            size_t h = cuts.at(edges["ac"]);  // ac的切割点
            
            // 检查切割点是否共面
            std::vector<Point3D> points = {nodes_[e], nodes_[f], nodes_[g], nodes_[h]};
            if (!checkCoplanar(points)) {
                return {tet};
            }
            
            return {
                // 三棱柱hgc-efd: acb-d方式切分
                Tetrahedron{h, c, g, e},
                Tetrahedron{e, f, c, g},
                Tetrahedron{e, f, d, c},
                // 三棱柱ahe-bgf: acb-d方式切分
                Tetrahedron{a, e, h, b},
                Tetrahedron{b, g, e, h},
                Tetrahedron{b, g, f, e}
            };
        }
        
        // 情况3：ab,bc,cd,ad被切割
        else if (cutEdges == std::unordered_set<std::string>{"ab","bc","cd","ad"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            size_t f = cuts.at(edges["bc"]);  // bc的切割点
            size_t g = cuts.at(edges["cd"]);  // cd的切割点
            size_t h = cuts.at(edges["ad"]);  // ad的切割点
            
            // 检查切割点是否共面
            std::vector<Point3D> points = {nodes_[e], nodes_[f], nodes_[g], nodes_[h]};
            if (!checkCoplanar(points)) {
                return {tet};
            }
            
            return {
                // 三棱柱cfg-aeh: acb-d方式切分
                Tetrahedron{c, g, f, a},
                Tetrahedron{a, e, g, f},
                Tetrahedron{a, e, h, g},
                // 三棱柱hgd-efb: acb-d方式切分
                Tetrahedron{h, d, g, e},
                Tetrahedron{e, f, b, g},
                Tetrahedron{e, b, d, g}
            };
        }
        
        else {
            return {tet};
        }
    }
    catch (const std::exception& e) {
        return {tet};
    }
}

std::vector<Tetrahedron> MeshCutter::splitCaseB(
    const Tetrahedron& tet,
    const std::unordered_map<Edge, size_t>& cuts) const {
    
    try {
        // 获取原始顶点
        size_t b = tet.v1, c = tet.v2, d = tet.v3, a = tet.v4;
        
        // 获取所有可能的边
        std::unordered_map<std::string, Edge> edges = {
            {"ab", Edge(a, b)}, {"ac", Edge(a, c)}, {"ad", Edge(a, d)},
            {"bc", Edge(b, c)}, {"bd", Edge(b, d)}, {"cd", Edge(c, d)}
        };
        
        // 检查哪些边被切割
        std::unordered_set<std::string> cutEdges;
        for (const auto& [name, edge] : edges) {
            if (cuts.count(edge)) {
                cutEdges.insert(name);
            }
        }
        
        // 情况1：ab,ac,ad被切割
        if (cutEdges == std::unordered_set<std::string>{"ab","ac","ad"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            size_t f = cuts.at(edges["ac"]);  // ac的切割点
            size_t g = cuts.at(edges["ad"]);  // ad的切割点
            
            return {
                // 四面体efg-a
                Tetrahedron{e, f, g, a},
                // 三棱柱efg-bcd: acb-d方式切分
                Tetrahedron{e, g, f, b},
                Tetrahedron{b, c, g, f},
                Tetrahedron{b, c, d, g}
            };
        }
        
        // 情况2：ad,cd,bd被切割
        else if (cutEdges == std::unordered_set<std::string>{"ad","cd","bd"}) {
            size_t e = cuts.at(edges["ad"]);  // ad的切割点
            size_t f = cuts.at(edges["cd"]);  // cd的切割点
            size_t g = cuts.at(edges["bd"]);  // bd的切割点
            
            return {
                // 四面体efg-d
                Tetrahedron{e, f, g, d},
                // 三棱柱efg-acb: acb-d方式切分
                Tetrahedron{e, g, f, a},
                Tetrahedron{a, c, g, f},
                Tetrahedron{a, c, b, g}
            };
        }
        
        // 情况3：ab,db,cb被切割
        else if (cutEdges == std::unordered_set<std::string>{"ab","bd","bc"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            size_t f = cuts.at(edges["bd"]);  // bd的切割点
            size_t g = cuts.at(edges["bc"]);  // bc的切割点
            
            return {
                // 四面体efg-b
                Tetrahedron{e, f, g, b},
                // 三棱柱efg-adc: acb-d方式切分
                Tetrahedron{e, g, f, a},
                Tetrahedron{a, d, g, f},
                Tetrahedron{a, d, c, g}
            };
        }
        
        // 情况4：ac,bc,dc被切割
        else if (cutEdges == std::unordered_set<std::string>{"ac","bc","cd"}) {
            size_t e = cuts.at(edges["ac"]);  // ac的切割点
            size_t f = cuts.at(edges["bc"]);  // bc的切割点
            size_t g = cuts.at(edges["cd"]);  // dc的切割点
            
            return {
                // 四面体efg-c
                Tetrahedron{e, f, g, c},
                // 三棱柱efg-abd: acb-d方式切分
                Tetrahedron{e, g, f, a},
                Tetrahedron{a, b, g, f},
                Tetrahedron{a, b, d, g}
            };
        }
        
        else {
            return {tet};
        }
    }
    catch (const std::exception& e) {
        return {tet};
    }
}

std::vector<Tetrahedron> MeshCutter::splitCaseC(
    const Tetrahedron& tet,
    const std::unordered_map<Edge, size_t>& cuts) const {
    
    try {
        // 获取原始顶点
        size_t b = tet.v1, c = tet.v2, d = tet.v3, a = tet.v4;
        
        // 获取所有可能的边
        std::unordered_map<std::string, Edge> edges = {
            {"ab", Edge(a, b)}, {"ac", Edge(a, c)}, {"ad", Edge(a, d)},
            {"bc", Edge(b, c)}, {"bd", Edge(b, d)}, {"cd", Edge(c, d)}
        };
        
        // 检查哪些边被切割
        std::unordered_set<std::string> cutEdges;
        for (const auto& [name, edge] : edges) {
            if (cuts.count(edge)) {
                cutEdges.insert(name);
            }
        }
        
        // 情况1：ab,ac被切割
        if (cutEdges == std::unordered_set<std::string>{"ab","ac"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            size_t f = cuts.at(edges["ac"]);  // ac的切割点
            
            return {
                // 四面体efd-a
                Tetrahedron{e, f, d, a},
                // 四棱锥befc-d切分为两个四面体
                Tetrahedron{b, e, f, d},
                Tetrahedron{b, f, c, d}
            };
        }
        
        // 情况2：bc,cd被切割
        else if (cutEdges == std::unordered_set<std::string>{"bc","cd"}) {
            size_t e = cuts.at(edges["bc"]);  // bc的切割点
            size_t f = cuts.at(edges["cd"]);  // cd的切割点
            
            return {
                // 四面体aef-c
                Tetrahedron{a, e, f, c},
                // 四棱锥befd-a切分为两个四面体
                Tetrahedron{b, e, f, a},
                Tetrahedron{b, f, d, a}
            };
        }
        
        // 情况3：ac,ad被切割
        else if (cutEdges == std::unordered_set<std::string>{"ac","ad"}) {
            size_t e = cuts.at(edges["ac"]);  // ac的切割点
            size_t f = cuts.at(edges["ad"]);  // ad的切割点
            
            return {
                // 四面体bef-a
                Tetrahedron{b, e, f, a},
                // 四棱锥cefd-b切分为两个四面体
                Tetrahedron{c, e, f, b},
                Tetrahedron{c, f, d, b}
            };
        }
        
        // 情况4：ab,ad被切割
        else if (cutEdges == std::unordered_set<std::string>{"ab","ad"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            size_t f = cuts.at(edges["ad"]);  // ad的切割点
            
            return {
                // 四面体ecf-a
                Tetrahedron{e, c, f, a},
                // 四棱锥ebdf-c切分为两个四面体
                Tetrahedron{e, b, d, c},
                Tetrahedron{e, d, f, c}
            };
        }
        
        else {
            return {tet};
        }
    }
    catch (const std::exception& e) {
        return {tet};
    }
}

std::vector<Tetrahedron> MeshCutter::splitCaseD(
    const Tetrahedron& tet,
    const std::unordered_map<Edge, size_t>& cuts) const {
    
    try {
        // 获取原始顶点
        size_t b = tet.v1, c = tet.v2, d = tet.v3, a = tet.v4;
        
        // 获取所有可能的边
        std::unordered_map<std::string, Edge> edges = {
            {"ab", Edge(a, b)}, {"ac", Edge(a, c)}, {"ad", Edge(a, d)},
            {"bc", Edge(b, c)}, {"bd", Edge(b, d)}, {"cd", Edge(c, d)}
        };
        
        // 检查哪些边被切割
        std::unordered_set<std::string> cutEdges;
        for (const auto& [name, edge] : edges) {
            if (cuts.count(edge)) {
                cutEdges.insert(name);
            }
        }
        
        // 情况1：ab边被切割
        if (cutEdges == std::unordered_set<std::string>{"ab"}) {
            size_t e = cuts.at(edges["ab"]);  // ab的切割点
            return {
                Tetrahedron{e, c, d, a},  // ecd-a
                Tetrahedron{e, d, c, b}   // edc-b
            };
        }
        
        // 情况2：cd边被切割
        else if (cutEdges == std::unordered_set<std::string>{"cd"}) {
            size_t e = cuts.at(edges["cd"]);  // cd的切割点
            return {
                Tetrahedron{e, a, b, c},  // eab-c
                Tetrahedron{e, b, a, d}   // eba-d
            };
        }
        
        // 情况3：bd边被切割
        else if (cutEdges == std::unordered_set<std::string>{"bd"}) {
            size_t e = cuts.at(edges["bd"]);  // bd的切割点
            return {
                Tetrahedron{e, c, a, b},  // eca-b
                Tetrahedron{e, a, c, d}   // eac-d
            };
        }
        
        // 情况4：bc边被切割
        else if (cutEdges == std::unordered_set<std::string>{"bc"}) {
            size_t e = cuts.at(edges["bc"]);  // bc的切割点
            return {
                Tetrahedron{e, d, a, c},  // eda-c
                Tetrahedron{e, a, d, b}   // ead-b
            };
        }
        
        // 情况5：ad边被切割
        else if (cutEdges == std::unordered_set<std::string>{"ad"}) {
            size_t e = cuts.at(edges["ad"]);  // ad的切割点
            return {
                Tetrahedron{e, b, c, a},  // ebc-a
                Tetrahedron{e, c, b, d}   // ecb-d
            };
        }
        
        // 情况6：ac边被切割
        else if (cutEdges == std::unordered_set<std::string>{"ac"}) {
            size_t e = cuts.at(edges["ac"]);  // ac的切割点
            return {
                Tetrahedron{e, b, d, c},  // ebd-c
                Tetrahedron{e, d, b, a}   // edb-a
            };
        }
        
        else {
            return {tet};
        }
    }
    catch (const std::exception& e) {
        return {tet};
    }
}

void MeshCutter::buildEdges() {
    edges_.clear();
    edges_.reserve(elements_.size() * 4);  // 使用更保守的预估值
    
    // 分批处理四面体
    const size_t batchSize = 1000;
    for (size_t i = 0; i < elements_.size(); i += batchSize) {
        size_t end = std::min(i + batchSize, elements_.size());
        for (size_t j = i; j < end; ++j) {
            const auto& tet = elements_[j];
            edges_.insert(Edge(tet.v1, tet.v2));
            edges_.insert(Edge(tet.v1, tet.v3));
            edges_.insert(Edge(tet.v1, tet.v4));
            edges_.insert(Edge(tet.v2, tet.v3));
            edges_.insert(Edge(tet.v2, tet.v4));
            edges_.insert(Edge(tet.v3, tet.v4));
        }
    }
}

void MeshCutter::buildEdgeToTets() {
    edgeToTets_.clear();
    
    // 分批处理四面体
    const size_t batchSize = 1000;
    for (size_t i = 0; i < elements_.size(); i += batchSize) {
        size_t end = std::min(i + batchSize, elements_.size());
        for (size_t j = i; j < end; ++j) {
            const auto& tet = elements_[j];
            std::array<Edge, 6> tetEdges = {{
                Edge(tet.v1, tet.v2), Edge(tet.v1, tet.v3),
                Edge(tet.v1, tet.v4), Edge(tet.v2, tet.v3),
                Edge(tet.v2, tet.v4), Edge(tet.v3, tet.v4)
            }};
            
            for (const auto& edge : tetEdges) {
                auto& tets = edgeToTets_[edge];
                if (tets.capacity() == 0) {
                    tets.reserve(4);  // 每条边通常属于2-4个四面体
                }
                tets.push_back(j);
            }
        }
    }
}

bool MeshCutter::checkCoplanar(const std::vector<Point3D>& points) const {
    if (points.size() < 4) return true;
    
    // 将Point3D转换为Eigen向量
    Eigen::MatrixXd coords(points.size(), 3);
    for (size_t i = 0; i < points.size(); ++i) {
        coords(i, 0) = points[i].x;
        coords(i, 1) = points[i].y;
        coords(i, 2) = points[i].z;
    }
    
    // 取前三个点计算平面法向量
    Eigen::Vector3d v1 = coords.row(1) - coords.row(0);
    Eigen::Vector3d v2 = coords.row(2) - coords.row(0);
    Eigen::Vector3d normal = v1.cross(v2);
    
    if (normal.norm() < 1e-6) {
        return false;  // 前三个点共线
    }
    
    // 检查其他点是否在平面上
    double d = -normal.dot(coords.row(0));
    for (int i = 3; i < coords.rows(); ++i) {
        if (std::abs(normal.dot(coords.row(i)) + d) > 1e-6) {
            return false;
        }
    }
    
    return true;
}

bool MeshCutter::isCoplanar(const std::vector<Point3D>& points) const {
    if (points.size() < 4) return true;
    
    // 将Point3D转换为Eigen向量
    Eigen::MatrixXd coords(points.size(), 3);
    for (size_t i = 0; i < points.size(); ++i) {
        coords(i, 0) = points[i].x;
        coords(i, 1) = points[i].y;
        coords(i, 2) = points[i].z;
    }
    
    // 计算中心点
    Eigen::Vector3d centroid = coords.colwise().mean();
    coords.rowwise() -= centroid.transpose();
    
    // SVD分解
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(coords, Eigen::ComputeFullV);
    
    // 计算奇异值比例
    const auto& s = svd.singularValues();
    double ratio = s[2] / s[0];  // 最小奇异值与最大奇异值的比例
    
    // 如果比例小于阈值，认为点集共面
    const double threshold = 1e-6;
    return ratio < threshold;
}

void MeshCutter::updateMeshData() {
    buildEdges();
    buildEdgeToTets();
}

bool MeshCutter::isEdgeInBounds(
    const Point3D& p1, 
    const Point3D& p2,
    const Point3D& boundsMin,
    const Point3D& boundsMax) const {
    
    // 检查边的包围盒是否与界面包围盒相交
    double minX = std::min(p1.x, p2.x);
    double maxX = std::max(p1.x, p2.x);
    double minY = std::min(p1.y, p2.y);
    double maxY = std::max(p1.y, p2.y);
    double minZ = std::min(p1.z, p2.z);
    double maxZ = std::max(p1.z, p2.z);
    
    return !(maxX < boundsMin.x || minX > boundsMax.x ||
             maxY < boundsMin.y || minY > boundsMax.y ||
             maxZ < boundsMin.z || minZ > boundsMax.z);
}

void MeshCutter::preprocessStrataInterface(const ImplicitSurface& surface, int surfaceType) {
    // 只处理地层界面，且跳过第一个(11)和最后一个界面
    if (surfaceType < 11 || surfaceType == 11 || surfaceType == lastStrataType_) {
        return;
    }
    
    // 获取界面的包围盒
    const auto& bounds = surface.getBounds();
    
    // 创建节点到四面体的映射，用于检查节点移动的有效性
    std::vector<std::vector<size_t>> nodeToTets(nodes_.size());
    for (size_t i = 0; i < elements_.size(); ++i) {
        const auto& tet = elements_[i];
        nodeToTets[tet.v1].push_back(i);
        nodeToTets[tet.v2].push_back(i);
        nodeToTets[tet.v3].push_back(i);
        nodeToTets[tet.v4].push_back(i);
    }
    
    // 遍历所有边
    for (const auto& tet : elements_) {
        std::vector<std::pair<size_t, size_t>> edges = {
            {tet.v1, tet.v2}, {tet.v1, tet.v3}, {tet.v1, tet.v4},
            {tet.v2, tet.v3}, {tet.v2, tet.v4}, {tet.v3, tet.v4}
        };
        
        for (auto& edge : edges) {
            if (edge.first > edge.second) {
                std::swap(edge.first, edge.second);
            }
            
            Point3D& p1 = nodes_[edge.first];
            Point3D& p2 = nodes_[edge.second];
            
            // 基本检查
            if (!isVerticalEdge(p1, p2) || 
                !canProcessEdge(p1, p2, surfaceType) ||
                !isEdgeInBounds(p1, p2, bounds.first, bounds.second)) {
                continue;
            }
            
            // 计算与界面的交点
            double phi1 = surface.evaluatePhi(p1);
            double phi2 = surface.evaluatePhi(p2);
            
            if (std::abs(phi1) < 1e-10 || std::abs(phi2) < 1e-10 || phi1 * phi2 >= 0) {
                continue;
            }
            
            // 计算交点
            Point3D intersection = computeCutPoint(p1, p2, surface);
            intersection.x = p1.x;
            intersection.y = p1.y;
            
            // 确定要移动的节点
            size_t nodeToMove;
            if (p1.type != 0 && p2.type == 0) {
                nodeToMove = edge.second;
            } else if (p1.type == 0 && p2.type != 0) {
                nodeToMove = edge.first;
            } else if (p1.type == 0 && p2.type == 0) {
                double dist1 = std::abs(intersection.z - p1.z);
                double dist2 = std::abs(intersection.z - p2.z);
                nodeToMove = (dist1 <= dist2) ? edge.first : edge.second;
            } else {
                continue;
            }
            
            // 检查移动节点是否会导致四面体翻转
            bool canMove = true;
            double oldZ = nodes_[nodeToMove].z;
            nodes_[nodeToMove].z = intersection.z;
            
            for (size_t tetIdx : nodeToTets[nodeToMove]) {
                const auto& tet = elements_[tetIdx];
                if (!isValidTetrahedron(tet)) {
                    canMove = false;
                    break;
                }
            }
            
            if (canMove) {
                nodes_[nodeToMove].type = surfaceType;
            } else {
                nodes_[nodeToMove].z = oldZ;  // 恢复原始位置
            }
        }
    }
}

bool MeshCutter::isVerticalEdge(const Point3D& p1, const Point3D& p2) const {
    // 检查两点的xy坐标是否相同（考虑浮点误差）
    const double epsilon = 1e-10;
    return std::abs(p1.x - p2.x) < epsilon && std::abs(p1.y - p2.y) < epsilon;
}

bool MeshCutter::canProcessEdge(const Point3D& p1, const Point3D& p2, int surfaceType) const {
    // 如果两个端点都是界面点，不处理
    if (p1.type != 0 && p2.type != 0) {
        return false;
    }
    
    // 如果任一端点是当前正在处理的界面点，不处理
    if (p1.type == surfaceType || p2.type == surfaceType) {
        return false;
    }
    
    return true;
}

bool MeshCutter::isValidTetrahedron(const Tetrahedron& tet) const {
    const Point3D& p1 = nodes_[tet.v1];
    const Point3D& p2 = nodes_[tet.v2];
    const Point3D& p3 = nodes_[tet.v3];
    const Point3D& p4 = nodes_[tet.v4];
    
    // 计算四面体体积
    double v321 = p3.x * p2.y * p1.z;
    double v231 = p2.x * p3.y * p1.z;
    double v312 = p3.x * p1.y * p2.z;
    double v132 = p1.x * p3.y * p2.z;
    double v213 = p2.x * p1.y * p3.z;
    double v123 = p1.x * p2.y * p3.z;
    
    double v421 = p4.x * p2.y * p1.z;
    double v241 = p2.x * p4.y * p1.z;
    double v412 = p4.x * p1.y * p2.z;
    double v142 = p1.x * p4.y * p2.z;
    double v214 = p2.x * p1.y * p4.z;
    double v124 = p1.x * p2.y * p4.z;
    
    double v431 = p4.x * p3.y * p1.z;
    double v341 = p3.x * p4.y * p1.z;
    double v413 = p4.x * p1.y * p3.z;
    double v143 = p1.x * p4.y * p3.z;
    double v314 = p3.x * p1.y * p4.z;
    double v134 = p1.x * p3.y * p4.z;
    
    double v432 = p4.x * p3.y * p2.z;
    double v342 = p3.x * p4.y * p2.z;
    double v423 = p4.x * p2.y * p3.z;
    double v243 = p2.x * p4.y * p3.z;
    double v324 = p3.x * p2.y * p4.z;
    double v234 = p2.x * p3.y * p4.z;
    
    double volume = (-v321 + v231 + v312 - v132 - v213 + v123
                    + v421 - v241 - v412 + v142 + v214 - v124
                    - v431 + v341 + v413 - v143 - v314 + v134
                    + v432 - v342 - v423 + v243 + v324 - v234) / 6.0;
    
    return volume > 1e-10;  // 使用一个小的正数作为阈值
} 